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COMPARISON OF A DISCRETE STEEPEST ASCENT METHOD 
WITH THE CONTINUOUS STEEPEST ASCENT METHOD 
FOR OPTIMAL PROGRAMING 

By A. Gary Childs 
Langley Research Center 

SUMMARY 

A discrete steepest ascent method which allows controls which are not piecewise 
constant (for example, it allows all continuous piecewise linear controls) is derived for 
the solution of optimal programing problems. This method is based on the continuous 
steepest ascent method of Bryson and Denham and new concepts introduced by Kelley and 
Denham in their development of "compatible" adjoints for taking into account the effects 
of numerical integration. The method is a generalization of the algorithm suggested by 
Canon, Cullum, and Polak with the details of the gradient computation given. The discrete 
method is compared with the continuous method for an aerodynamics problem for which 
an analytic solution is given by Pontryagin's maximum principle, and numerical results 
are presented. The discrete method converges more rapidly than the continuous method 
at first, but then for some undetermined reason, loses its exponential convergence rate. 

A comparison is also made for the algorithm of Canon, Cullum, and Polak using piece - 
wise constant controls. This algorithm is very competitive with the continuous algorithm. 

INTRODUCTION 

One of the simplest and most direct algorithms for the solution of optimal programing 
problems is the steepest ascent algorithm of Bryson and Denham. (See ref. 1.) This 
method, being a first-order method, suffers from poor terminal convergence. Kelley and 
Denham (ref. 2) claim convergence can be improved for a conjugate gradient algorithm by 
use of "compatible" adjoints. These adjoints satisfy adjoint difference equations for the 
numerical integration difference equations associated with the differential constraints and 
represent an attempt to take into account numerical calculations in the algorithm. In this 
paper, what is called "a discrete steepest ascent" algorithm is derived; this algorithm is 
based on the numerical integration difference equations and a reformulation of the optimal 
programing problem as a discrete problem. This algorithm is a generalization of the 
algorithm given by Canon, Cullum, and Polak in reference 3. The new algorithm is then 
compared with the Bryson-Denham algorithm for an aerodynamics problem for which 


an analytic solution is determined by Pontryagin's maximum principle. (See ref. 4.) 
Numerical results comparing the convergence rates of the two methods are presented 
in graphical form. Also, the algorithm of Canon, Cullum, and Polak is compared with 
the continuous algorithm for piecewise constant controls. Numerical results comparing 
convergence rates are presented. 


SYMBOLS 

Measurements and calculations were made in U.S. Customary Units. They are 
presented herein in the International System of Units (SI). 


C D 


drag coefficient 


CpQ L 2 induced drag coefficient 

Cd,o zero-lift drag coefficient 

C L lift coefficient 


C L 


a 


D 


lift -curve slope 


drag, newtons 


d/3 = drf/ - 

dS2, defy, di// defined by equations (7) 
N* 

(dP) 2 = ^ Sufw^uj 
i=0 


& 


function of Xj,Uj,Uj. i,tj,t^ + i equal to Axj 


T'x - '^( x i> u i> u i+l>* : i>t : i+l) 


Fi = 




1 8xj 


f _ 9f 
* x " ax 


fimetion of x(t), u(t), and t equal to x(t) 
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defined by equations (Cl) or (C2) 



P 0 ,P 1 ,P 2 Pontryagin adjoints 

n 

S aerodynamic surface area, meters* 

T thrust, newtons 

t time or independent variable, seconds 

tf final time, seconds 


3 



initial time, seconds 


u 


u r 


U N 


control vector or horizontal velocity, meter s/second 


initial control at initial increment 


final control at final increment 


V = (u2 + v2) ^ , meters/second 


W 


Wi 


vertical velocity, meters/second 


weight, newtons 


weighting matrix (symmetric) 


state vector 


initial state 


a = 9 - y, radians 


r 


flight -path angle, radians 


control angle, radians 


difference -equation adjoint vector 


x Q,,i’ X \p, i> 

X <M 


difference -equation adjoints associated with O, c 


x 4>J2,i> x ^S2,i defined by equations (9) 
p, v Lagrange multipliers 

p air density, kilograms/ meter^ 


4 > 


performance index for development of algorithm 


i, and i// 
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ij; constraint vector function 

SI stopping condition function 

Operations: 

6( )i = ( )i - ( )j* 

A( H-Ow-Oi 


4 

definition 


• 

( ) 

derivative of ( ) with respect to 

t 

d( ) 

differential of ( ) with respect to 

^ 6u o\ 



^ 6u N/ 

n - 1 

matrix inverse 



!□ 
9 ( ) 

() T 


partial derivative of 3 with respect to ( ) 
matrix-vector transpose 


A symbol with an asterisk denotes a nominal variable; a symbol without an asterisk 
denotes a general or arbitrary variable. 

ANALYSIS 

The Optimal Programing Problem 
Let the following problem be considered. Maximize 


^(tf^tf) 


subject to 


x = f(x,u,t) 


(to St St,) 
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and 

= 0 

with tf determined by 

= o 

and t 0 and x 0 as given constants. For this problem, <fi is a scalar function called 
the cost function, £2 is a scalar function called the stopping condition, x is a vector 
function called the state, u is a vector function called the control, and i// is a vector 
function called the terminal constraint. The dimensions of the vector functions are, in 
general, different. The optimal programing problem in this form is discussed in refer- 
ence 1 along with the continuous steepest ascent method for solving it. 

The theory of reference 1 requires (strictly speaking) that all functions be known 
exactly. However, for the solution of the foregoing problem, one frequently uses numer- 
ical integrations for obtaining the state. Hence, an algorithm is required which recog- 
nizes the use of difference equations for the state, that is, a discrete algorithm. A dis- 
crete method is developed in an analogous fashion to the development of the continuous 
method. A reformulation of the problem taking into account some of the before -mentioned 
realities of computer calculations is made. Maximize 

0( x n^n) 

subject to 


Axj — Xj^j - Xj = 2^(xj>Uj,Uj^i,t^,tj^j^ A (i — 0, . . ., N - 1) (1) 

and 

^( x n^n) = 0 


with tjj determined by 

^(xn^n) = 0 

and t Q and x 0 as given constants. For this problem, <p, £2, Xj, . Uj, and \jy are 
scalars or vectors as in the continuous formulation of the problem. In addition, it is 
specified that 

ti+l - ti = h (i = 0, . . ., N - 2) 
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where h is a constant. Also, 


*N " %-l = h 

Hence, the t^ are the time points of the numerical integration of the continuous state 
equations (represented by eq. (1)) with tjq- = tf. 


Formulation of Discrete Method 

Let uj* (i = 0, . . ., N*) be a nominal control. Let x * (i = 0, . . ., N*) be the 
nominal state resulting from this control. Suppose a new control Uj (i = 0, . . ., N*) is 
given which yields a new state (i = 0, . . ., N*). Then 


6uj A Uj[ - uf 
6x i = x i “ x * 


(i = 0, . . ., N*) 


Expanding equations (1) about the nominal control and state and truncating to the linear 
terms yield 

A(6xj) = 6x i+ j - 6xj ® FjSxj + G^6u^ + HjSu^j (i = 0, . . ., N* - 1) (2) 


where 


G i A |^ x l» u l» u i+l» t i; t i+l) 
H i £ Uf, Uf + 1, ti,t i+ l) 




The adjoint equations for equations (2) are obtained from 

A ( x ? 6x i) = x l 


v i+l 6x i+l " X ? 6x i 


= x ?+l6x i+ i - xE-lSxi + xf + i6Xi - xffiXi 

= X^jA^Xj) + A^X^6x^ 

x {+l F i + A ( x i’)] 6x i + X f+l( G i 6u i + H i 6u i+l) 


(i = 0, . . ., N* - 1) 


(i = 0, 


., N* - 1) 


(3) 
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Then defining the Xj terms to be the compatible adjoints of reference 2 for the partic- 
ular type of numerical integration used in equations (1), 


A(A|) — — -Fj Xj^j (i — 0, • . N 1) 

and substituting equations (4) into equations (3) gives 

A(x? , 8x i ) * X^jfGiSUi + H^u^) 

Summing equations (5) from i * 0 to N* yields 

N*-l 

Xjq-^Sx^* ~ Xq8x q + ^ Xj + j(G^6uj + H^6U| + j) 
i=0 

Linear approximations for changes of \f/, and f2 are 

d< ^ “ ^( X N* ,t N*) 6x N* + ai4i( X N* ,t N*) f ( X N* ,U N* ,t N*) + %f( X N* ,t: N*) (*f " 4 N*) 


( 4 ) 


(5) 


( 6 ) 


a<p 

dx 


^N^N*) 6x N * + -J±L(x N *,t^ f(x N *,u N *,t N ^ + ^£-(x N *,t N *j ^t f - t N *j 
^( X N*’ t N*)( X ( t f) " X N*) + at^( X N* 5t N*) (*f " *N*) + ^( x N*» t N*) ( X N* ~ X N*) 


ax 


N' 


~ 0(x(tf),t f ) - 0(x N *,t N+ ) + 0(x N *,t N *) - cf)(x**,t N *) 

= 0(x(tf),tf) - ^(x**,t N *) 

— — 
- al^N^N*) 6 x N* + 9^;( x N*» t N*) f ( x N*’ U N*’ t N*) H "^( x N*’ t N*) 

* ^( x {hyt) ~ ^( x N* ,t N*) 


dft = ^^( x N* ,t N*) 5x N* + |axJ J ( X N*’ t N*) f ( x N*’ u N* ,t N*) + ^( x N*’ t N*) (*1 " l N*) 


\ (?) 


(t f - t N *) 
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where tf is the final time associated with the control uj (i = 0, . . N*) and is deter- 

mined by use of the stopping condition S2. If 0 is not zero by the time t-j^, the new 
control must be extrapolated in some way and the numerical integration of the state con- 
tinued until a zero of O is found. Note that a reversion to the continuous system 
(x = f(x,u,t)) is used to obtain an approximation for x^ - x-^ *. 

Let 


X 4>,N* = '^( x N*’ t N*) 

X J,N* - ^( x N* t N*) 

X «,N % - ^( X N^*N*) 
and using equation (6), equations (7) become 


N*-l 




d 4> = X <£,0 6x o + £ x cj),i+l( G i 6u i + H i6u i+1 ) 


i=0 


3 <p 
3x 


( x N +,t N*) f ( x N*’ u N*» t N*) + ( 4 f ' *N*) 


N*-l 


d^ = xJ, 0 5x o + J xJ >i+ i(G i 6u i + Hi5u i+1 ) 


i=0 


^( x N*’ t N*) f ( x N* ,u N* ,t N*) + 3^( X N*» t N*) ( fc f " *N*) 


N*-l 


dS2 = X n,0 6x o + X X ^,i+l( G i 6u i + H i 6u i+l) 


i=0 


_3S2 

3x 


^( x N*’ t N*) f ( x N*» u N*> t N*) + l^( X N* ,t N*) ( t f " *N*) 


( 8 ) 
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Setting dfi = 0 gives a formula for tf - t N * which modifies equations (8) to 


N*-l 


d< £ = ^OjO^o + ^ A 0J2,i+l( G i$ u i + **i^ u i+l) 


i=0 


N*-l 


T <p 

d ’/ / = V«,0 6x o + 2j X W+l( G i 5u i +H i 6u i+l) 


i=0 


0) 


where 


.T a\T 

= A 0,i - m 


f ( x N*’ u N*» t N*) + 9t^ : ( X N*’ t N*) 
^~( x N*» t N*) f ( x N 5 ^ u N*» t N*) + ^( X N*» t N*) 


l S2,i 


> (i=0, 


,T A> t 
A \pC2,i = A i//,i “ 


— f ( x N^ u N*>t N *) + it^( X N* ;t N*) 

^( X N**N*) f ( x N*’ u N*» t N*) + %( X N* t N*) 


_an 

8x- 


S2,i 


N*) 


It is desired to maximize 




T 

<p£l,i+l 


^ Tx Jfi,i + l) G i 


+ ( x 0J2,i " yTx Jfi,i) H i-l " MSuJ’Wi 


6uj + v T dip + ^(dP)' 


with respect to Su 4 (i = 0, . . N*) where 6x 0 , dip, and the step size 


N* 

(dP) 2 4 ^ euJ’WjSuj 
i=0 


(10) 


(where is a symmetric m x m matrix) are assumed to be specified and v and p 
are Lagrange multipliers. Also, Cr N * and H_j are defined to be zero matrices in con- 
structing a form for d <p which will be easier to maximize. Taking the differential yields 
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N r \ 

d(d 0) = £ ( x 0S2,i+l - ^Jo,i+l) G i + ( x 0f2,i - ^Jn,i) H i-l - 2/iSufWi 


d (6Ui) 


This relation is zero for all values of (i = 0, . . N*) if 

5u i * sir w r 1 [ G i : ( x * o , i + i - + * f - i (\ fo,i - vn , i -) 


(i = 0, . . N*) (11) 


Substituting these control increments into the equations for dp and (dP) 2 gives 
values for v and ji. By use of these values, the 6uj (i = 0, . . N*) are completely 

specified. This process is carried out in appendix A. The result is 


6uj = 


-|l/2 

(dP) 2 - d/ 3 T I^ d/3 1 w _ j 


*00 I 00 I 00 I 00 ^ 


(I + + iq.! 


lT 


X 0£2,i " *-p£l,i^pp*-pp 


+ Wj 


(i + F t ) 1 G i + H i _ 1 


-iT 


X pto,i i: pp 


(i = 0, . . N*) (12) 


where 


N T 

V* = I ( l + f i)' 1g i + Hi-i w f 1 (i + F i)" lo i + H i-i| > 


i=0 

N* 




i=0 
N* 


r nr ~iT 

V<#> = Z ^O.i? + F i)" 1(3 i + Hi-iki'p + f i)' 1g i - H i-iJ x w ) 


n r -i r -rT 

J 00 = Y x 0fi,i i 1 + F i)’ lG i + H i-1 W i 1 (* + F i)’ lo i + H i-1 > 


i=0 


l 0O,i 


(13) 


The change of 0 caused by this control change is (by use of eq. (4)) 

T * 


N 

d 0 = X 0S2, 0 6x o + Y X 0«,i l 1 + F i) G i + H i-1 
i=0 


6u, 


"1 


1/2 


= X 0S2,O6^O + |l_( d P> 2 - d ^)(l00 - ij + 'Ifyp ^ 


(14) 
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The Sample Problem 


A problem was chosen for the comparison of the usual steepest ascent method with 
the discrete steepest ascent method. The final vertical velocity of an airplane with aero- 
dynamic parameters similar to those of the X-15 was to be maximized. The airplane 
was flying in a vertical plane, that is, it was constrained to two-dimensional motion. It 
had no initial horizontal velocity and some positive vertical velocity. The airplane was 
considered to be a point mass. The coordinate system and force diagram used for the 
problem are shown in figure 1. The equations of motion are 

mu = T cos 9 - D cos y - L sin y 


mv = T sin 9 - W - D sin y + L cos y 


where T is the (constant) thrust, W = mg is the weight, D is the drag, L is the 
lift, m is the mass, g is the earth's gravitational acceleration constant, u is the 
horizontal speed, v is the vertical speed, 9 is the control angle, and y is the flight- 
path angle. The formulas used for the lift and drag are 


L = C L 


D = C d 


pV 2 S 

2 

pV 2 S 

2 


where 


C 


L 


= C 


a 


a 


Cr» = Ct 


'D ~ ^D,o + C D C 2 (J l> a a ‘ 


j/2 

and p is the air density (constant), V = (u2 + v2) / is the speed, S is the wing sur- 
face area, and Cj^ , Cp Q , and Cj^ g are constant coefficients determined by the aii 

plane configuration. The values of constants used in the computer study are as follows: 


T, N . . . 
m, kg . . 
g, m/sec2 
Ct . . . 

■*— *Q£ 

c D,o • • * 

Cd C l 2 • 

p, kg/m 3 
S, m2 . . 


250 000 
15 000 
9.7759 
3.61 
0.055 
0.398 
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0.41351 

20 



The problem, then, was to maximize v(tf) - v(t Q ) where tf was given and u(t 0 ) 
was zero, v(t 0 ) = v 0 > 0. This problem was solved analytically by using Pontryagin’s 
maximum principle. (See ref. 4.) The solution was to fly with 9 constantly equal 
to 90°. This solution is derived in appendix B. The problem was programed for solution 
by both the usual steepest ascent method and the discrete steepest ascent method. Fourth- 
order Runge-Kutta integration was used (ref. 5), and appendix C gives the relationship 
between the ad joints of the two methods. Linear interpolation was used to obtain inter- 
mediate values of the control. A comparison of the results from the two methods obtained 
by the computer is given in the next section. 

RESULTS AND DISCUSSION 

Computer solutions of the problem were programed for both steepest ascent algo- 
rithms. The discrete algorithm compared very unfavorably with the continuous algorithm. 
The discrete algorithm gave corrections for 9 Q and which are only about half as 

much as the corrections given by the continuous algorithm. Recall that H_j and G^* 
in equation (12) are zero matrices. Hence, there are only half as many terms in the com- 
putation for 5u Q and 6u N * as for the other (i = 1, . . ., N* - 1). As can be seen 
in reference 1, for the continuous algorithm, the equation corresponding to equation (12) 
has an equal number of terms for the calculation of all 6u(t) (t Q S t S tf*). If one traces 
back both from equation (12) and the corresponding equation in reference 1, one discovers 
that the difference just noted in these equations appears first in the state equations. In 
this paper, the equations are 

AXi = F(x i ,u i ,u i+1 ,t i ,t i+1 ) (i = 0, . . ., N* - 1) 

In reference 1, the equations are 

x = f (x,u,t) (*0 = * = tf) 

The difference arises because both u^ and Uf + j appear in the first set of equations. 
Thus, in summing to get x N *, the uj (i = 1, . . ., N* - 1) appear twice as often as u Q 
and u N *. For the continuous equations of motion, every u(t) (t Q S t § tf) has equal influ- 
ence in determining x^tf). Since the first set of equations represents a numerical inte- 
gration process for the second set, the difficulty appears in the attempted discrete 
modeling of a continuous model for the dynamics of the problem. The most obvious solu- 
tion for such a difficulty is to find directly a discrete model for the dynamics of the prob- 
lem. This model would be a set of difference equations of the form 
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that is, equations (1) where it is assumed that uj + i and tj + ^ do not appear. 

This form of the difference equations for the state is essentially the same as that 
used in reference 3. It can be achieved by requiring the control of the sample problem to 
be constant over the numerical integration intervals. In general, however, there might 
exist "wild" piecewise constant controls (that is, where the control would have, at some 
interval end points, large jump discontinuities) which, although optimal, would be imprac- 
tical to implement. Also, even though such controls, as optimal controls, might not occur, 
the possibility exists that such controls would be encountered in the iterations of the algo- 
rithm and could precipitate divergence. Putting u^ + j and tj + j in the difference equa- 
tions allows continuity restrictions to be used to keep the control reasonable. The fore- 
going should be considered an argument justifying the study of the algorithm of this paper 
rather than as a criticism of the effectiveness of the algorithm of reference 3. Indeed, 
results presented for this algorithm in figure 2 reveal that it is very competitive with 
the continuous algorithm for this problem where the restriction has been made for both 
methods that 6 be piecewise constant. Appendix D gives the derivation of the equations 
used in this study which was suggested by Terry A. Straeter of the Langley Research 
Center. 


To attempt an alleviation of the before -mentioned difficulty, an empirical device 
was employed. The corrections 6u Q and 6u N * were doubled. It was hoped that this 
doubling would give 6u Q and 6u^* an almost equal influence with the other 8u^ terms. 
For the sample problem, 60 o and 60 N * would be doubled and would be about right. 

This procedure helped, but performance of the discrete algorithm still was not as 
good as that of the continuous algorithm. Further computer experience revealed that the 
6 terms (i = 0, . . ., N*) were larger for the discrete algorithm than for the continuous 
algorithm. The formula for these corrections for the sample problem is 


60 i = 


(dp y 


00 


1/2 j 


(i + F^Gi + Hj.j 


uT 


'•0J2,i 


Observe that an increase in I^ would make these quantities smaller. On the basis of 
the problems with 56 0 and 69 -$*, examination of the equation for 1^,0 of equations (13) 
revealed that there are only half as many terms for i = 0 and i = N* as there are for 
the other values of i. Therefore, the i = 0 and i = N* terms were doubled. After 
these two changes, the performance of the discrete algorithm modified in this fashion 
during the first few iterations, was better than that of the continuous algorithm. 
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Figures 3 to 9 give the results of a comparison of the two algorithms as they are 
presently formulated. Each case was run for 25 iterations. Points are omitted for iter- 
ations where the automatic convergence scheme rejected a forward trajectory which did 
not give as good a performance index as the trajectory of the previous iteration. 


For figure 3, the final time was 10 seconds, the computation interval was 0.5 sec- 
ond, and the nominal control was 0 s 70°. The computation interval of 0.5 second pro- 
duces terms twice as large for the discrete (dP) 2 as would appear in an approximation 
for the continuous (dP) 2 ; that is, 

N* 

Discrete: (dP) 2 4 Y 60? 


i=0 


Continuous: 



N*-l 

60 2 dt ~ ^ 
i=0 


1 

2 



Hence, the initial (dP) 2 for the discrete algorithm was chosen to be twice the value of 
the initial (dP) for the continuous algorithm to allow an equal amount of control change 
for both algorithms. 

For a fixed -final -time fixed-initial-condition problem, 


d0 = [l^(dP) 2 ] 1/2 

Therefore, an I ^ for the discrete algorithm equal to half the value of the for 

the continuous algorithm will give just as much change for cf>. Hence, the ordinates for 
figure 3 are I ^ for the continuous method and 21^ for the discrete method. Sim- 
ilar adjustments were made for the results shown in figures 4 to 9. 

The case of figure 3 was considered as the base problem. Figures 4 to 9 show 
results of perturbations from this case. These cases were run to determine whether any 
changes of certain parameters would improve the performance of the discrete algorithm 
relative to the continuous algorithm. In figures 4 and 5 results are presented for per- 
turbations in the computation interval to 1 second and 5 seconds, respectively. In fig- 
ures 6 and 7 results are given for perturbations from the base problem nominal control 
to 0 s 80° and 0 s 60°, respectively. The values of initial (dP) 2 for these two cases 
were chosen to give less and more control effort to effect convergence in about 25 itera- 
tions. Figures 8 and 9 give results for final times of 20 seconds and 30 seconds, respec- 
tively. For both cases the computation interval was 5 seconds. This calculation was an 
attempt to determine the effect of large errors in the numerical integration. 
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Both methods converge for the base problem. As shown in figure 3, the discrete 
method converges faster than the continuous method until it loses its exponential con- 
vergence rate for some unknown reason. The loss of convergence rate is seen as a 
"flattening out” on the semi -log plot. This behavior is typical of the discrete method 
and is seen in the other figures where the discrete method converges. No good hypothe- 
sis has been formulated to explain the behavior. It is noted, however, that the discrete 
method stops working in each case after about the same number of iterations. An exami- 
nation of the computer output reveals that the control after this number of iterations is 
"rough;" that is, the values of the % are scattered on both sides of the theoretically 
optimal 90°. For the continuous method, the control is kept relatively "smooth" in the 
same sense. Hence, the reason for the loss of convergence of the discrete algorithm may 
be involved with some implicit constraints which the continuous algorithm exerts on the 
control to keep it smooth. This reasoning leads one to think about the possibility of some 
uniqueness problem for the discrete algorithm (caused by the lack of the implicit con- 
straints), especially since such problems are commonly encountered in dealing with dis- 
crete models for continuous systems. Another possible explanation for the flattening out 
is that it is due to the modifications of the discrete algorithm described in the first part 
of the section "Results and Discussion;" that is, a switch back to the original algorithm 
when the flattening out begins may improve convergence. Unfortunately, time limitations 
prevented a study of this possibility. 

Comparison of figures 3 to 5 reveals the effect of a change in the computation inter- 
val. Both methods converge more slowly for a larger computation interval. The discrete 
method is notably more sensitive to a degrading of numerical accuracy in the integration, 
particularly in the flat section. A possible reason for this greater sensitivity is that as 
the final time remains the same and the computation interval increases, the number of 
integration points becomes less. Hence, the problems mentioned earlier with respect to 
50 o and have a greater influence on results. The computation interval of 5 sec- 

onds, yielding only three integration points, affected the discrete algorithm so adversely 
that it diverged. 

Comparison of figures 3, 6, and 7 demonstrates an interesting phenomenon. The 
flat part of the discrete algorithm plot occurs at different levels for different choices of 
the nominal control. It is especially surprising that the convergence "settles out" for a 
60° nominal at a lower value of the gradient than for a 70° nominal. The convergence 
rates with this exception are about the same for different nominal controls; that is, the 
average slope of a curve connecting the points would be about the same for the three plots. 

The discrete method diverged for figure 8. The continuous method was converging 
very slowly for the case of figure 8, as is indicated by the presence of points for up to 
24 iterations. Both methods diverged in figure 9. The erratic behavior of the points 
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represents an inaccurate calculation of the gradient and explains the poor convergence 
rate. The state integration up to 10 seconds is accurate to three significant figures for 
the case of figure 8. The same would, of course, be true for the case of figure 9 since 
the computation interval is the same as that for figure 8. However, the error over the 
30 seconds of figure 9 would be greater than the error over the 20 seconds of figure 8. 
Therefore, the state integration error probably contributed substantially to the divergence 
in these cases. 


CONCLUDING REMARKS 

A discrete steepest ascent algorithm which takes into account numerical integration 
of the differential constraints and allows controls which are not piecewise constant has 
been formulated for the solution of optimal programing problems. The algorithm was 
compared with the continuous steepest ascent algorithm of Bryson and Denham for an aero- 
dynamic problem for which an analytic solution had been obtained by using Pontryagin's 
maximum principle. For this problem the discrete algorithm converged somewhat faster 
initially but eventually slowed its convergence rate greatly. Basic difficulties with the 
discrete method which caused this condition are described. The prime and seemingly 
unavoidable difficulty is that the method of this paper uses a discrete model for a problem 
stated in terms of a continuous model. For physical problems formulated by use of such 
a continuous model, the method would apparently be of limited usefulness. However, if 
the problem is stated in terms of a discrete model (that is, using piecewise constant con- 
trols), the algorithm of Canon, Cullum, and Polak, of which the algorithm of this paper is 
a generalization, is very competitive with the continuous algorithm. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., October 13, 1971. 
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APPENDIX A 


DERIVATION OF EQUATIONS FOR CONSTANT 
LAGRANGE MULTIPLIERS 

In this appendix formulas are determined for the Lagrange multipliers, v and p. 
introduced in the section "Analysis." The control corrections 6uj are then given by 
using these formulas. Using equations (9) and (11) for di}/ and 6uj (i = 0,. . .,N*) 
yields 

N* 

w = x lto, 0 6x o + Z Wfi,i+l Gi + X W H i-l )^r W i _1 R ( X 0O,i+l " V.i+l 1 ') 


i=0 


+ H L-1 ( x 0O,i ■ ^0,1*')] ■ X i//fi,0 6x o + 2p (v<£ ‘ W*) 

where (by using eqs. (4)) 

N* 

V0 = Z ( X Js2,i+l G i + x Jo,i H i-l) W i 1 ( G ? X ^0,i+l + H i-l X 0fi,i) 

i=0 

- f + F i)" lG i + ^-l| w i ^ + F i)' lQ i + H i-l] T Vo,i 


(Al) 


N* 

W = Z ( X Jn,i+l G i + ^^^,i H i-l) W i 1 ( G i , ^^fi,i+l + H i-l^^n,i) 


i=0 

N 


- Z "kiji 1 + F i)“ lQ i + H i-i] w i "E 1 + F i)" lo i + H i-i] V«,i 


Solving for ^ from equation (Al) yields 
- - - 2(X <») 


where 


(A2) 


& = - X W,0 5x o 
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APPENDIX A - Continued 


Using equations (10) and (11) for (dP)2 and 6uj (i = 0, . . N*) and equation (A2) 


N* r 


( dp ) 2 = ~2 ^ |(^S2,i+l " v Tx $S2,i+l) G i + (*-$S2,i - v T $Si,i) H i-1 


' X xp^,i+l v ) + H ?-l(^0O,i - x ^,i v ) 

- 2vTl n + ‘' Tl «'") 


w 


-1 


^0S2,i 


i+1 


4 M 


4jL^ 


~<p<p ~ ^{^xp<p ~ 2 ^ ) ^xf/xjAxpcj) + ^ (^">p(f> ~ 2 ^ d ^) 


4, 2 (^0 0 " *4 / (p^4 / *''l /( P + ^ ^xf/xf/ 


XI 

*$(/> = ^ (^0f2,i+l^i + ^$S2,i®i-l)^i (®i X (p£l, i+1 + **i-l^</>fi,i) 


4j a 

where (by using eqs. (4)) 

N* 

I 

i=0 
N* 
i=0 

Solving for 2/z from equation (A3) yields 

4/2 


(A3) 


2 x w 


( \-l 

_ 1 

/ \-l 

(l + Fi) Gi + Hi.i 

Wi 1 

(i + Ff) Gi + Hi.j 


l 0n,i 


2/i = 


T -1 

I<P4> ~ 


(dP) 2 - d/3 T I^, d/3 


(A4) 


where the positive sign has been chosen to make dtp positive. 
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APPENDIX A - Concluded 


By using equations (11), (A2), and (A4), the control change equations are given by 
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APPENDIX B 


ANALYTICAL SOLUTION DERIVED BY USE 
OF PONTRYAGIN'S PRINCIPLE 

Pontryagin's maximum principle (ref. 4) is applied to the same problem introduced 
earlier. It is established that the control 9 = 90° satisfies a necessary condition for 
optimality. Let 


K D,1 ^ (f)c D ,o 

K D ,2 £ (£)c Dcl 2C l a l 


k t a 

L * 


fps'lcx 
\2 r h a 


(Bl) 


From figure 1, it is clear that 


cos y = ^ 




Sin y - 


a = 9- y= 6- cos"* ^ 


(0 S y S n) 


(B2) 


where a is the angle of attack of the airplane. By using equations (Bl) and (B2), the 
equations of motion can be rewritten as 


a/2 


= T cos 9 - -L + Kj) 2 a ^) u + Kl«vJ(u 2 + 
mv = T sin 9 - mg - jj^K D i + 1% 2 « 2 )v “ K L au|(u2 + v 2 )^ 2 


(B3) 


It is desired to minimize 


J = V - irdt = -[ v ( t f)- v (to)] 
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APPENDIX B — Continued 


The Pontryagin pseudo-Hamiltonian is 


H § p Q (-v) + pjii + p 2 v = Pg m P ° <|r sin 9 - mg - [(%,! + Kb,2 a2 ) v “ K L a uj( u2 + v2) 1//2 


^<4 cos 9 - [(%,! + k D,2« 2 ) u + K L av|(u2 + v 2) 1 / 2 | 


where 


Pi £ " Sf = ‘ mV | 2K D,2“ I f) u + ( Kd 4 + K D,2“ 2 ) + k L § v _ ( u2 + v2 ) 1/2 


- ( k D, 1 + K D,2 a2 ) u + K L avJu(u2 + v2)' 


■1/21 _ P 2 ~ p 0 / 1/W, oQ! 9a\ 

du) 


m 


(( 2k d,: 


' K L |j u " K l cJ(u 2 + v2)!/2 - [(k d>1 + K d 2 (» 2 )v - K L auJu(u2 + v 2 ) 
P 2 $ ' f = - Hrf[( 2K D,2“ If)" + k L §7 v + K L «]( U 2 + v2)l/2 


- ( k D,1 + k D,2“ 2 ) u + K L avJv( U 2 + v2)‘ 1 A - P2 m P ° |- (2% ,2“ ff)v 
+ ( K D,1 + K D,2 a2 ) " K L If u ]( u2 + v2 ) 1//2 - ( K D,1 + K D,2 “ 2 ) v 


-iM 


> (B4) 


- K-£au|v(u2 + v2) - ^/ 2 ^ 


By using the equation for a from equations (B2), 

da _ \Jy 2 
9 U "y 2 


3a = UV 

8v ‘ ’ pv 2 
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APPENDIX B - Concluded 


One has also 



Pift) ■ 0 


p 2 (*f) = 0 

From reference 4, in order to minimize J with respect to 0, H must be maxi- 
mized with respect to 9. Therefore, 

|f = ^{' T sin 0 - [( 2K D,2 «) u + %/|(u2 + vZ) 1 /^ + £2_£0 ^ CQS 0 

' [( 2K D,2“) v ‘ K L u ]( u2 + v 2) 1/2 J =0 (t Q § t £ tf ) 

The claim is made that 9 = |- (t Q = t ^ tf) will be a solution. By putting 9 = |- 
into the equation for u of equations (B3), it is seen that u = 0 ^t Q ^ t i tf) satisfies the 
equation. This result is in agreement with the boundary condition u(t Q ) = 0. From the 
equation for a from equations (B2), a = 0 ^t Q = t § tf). By putting 0, = |- and u = 0 
into the equation for p^ of equations (B4), Pj = 0 (t Q S t g tf ) is seen to be a solution 
which also satisfies p^tf) = 0. With 0 = £ u = 0, and p.^ = 0, it is seen that 

f-ofcit.n). 

The second partial is 

= 5f(- T cos 9 - [( 2k D, 2) u ]( u2 + v 2 ) 1/2 } + V*[-T sin B - [(2K D 2 )v](u2 + v2) l/2 j 
- ^(- T - 2 %,2^) 

where 9 - !■, u = 0, p n = -1, and p 2 = 0 have been used. The value will be less 

1 u 1 802 

than zero for t Q S t = tf if p£ + 1 > 0, and -T - 2 Kj- ) ^ v /v 2 will be less than zero 
for t Q S t 2 tf. These two conditions can be shown to hold if it is further assumed 
that T > mg + KpfV^ and that v is continuous. Thus, under these conditions, 

$ *= g- (t Q = t = tf) satisfies a sufficient condition for a local maximum of H and 
thereby satisfies a necessary condition for optimality. 
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APPENDIX C 


DERIVATION OF COMPATIBLE AD JOINTS FOR FOURTH- ORDER 
RUNGE-KUTTA INTEGRATION 

In this appendix the adjoints for both the continuous method and the discrete method 
for fourth-order Runge-Kutta integration are developed. Relationships between the two 
methods are discussed. For both methods x* is given by Runge-Kutta integration of 
the equation 

X* = f(x*,u*,t) 

with initial condition x*(t Q ) = x Q and the interval t 0 = t = t£ is partitioned as in the 
section "Analysis." 

The adjoints for the usual steepest ascent method are generated by 
X = -fJ(x*,u*,t)A 


with the terminal condition A(tf) for <fi and \p the same as for the discrete method. 
For any i = 0, 1, . . ., N* - 2, by using backwards Runge-Kutta fourth-order integration, 


M - H+l * j}fe<*>*Ul * *£<»(»i + l + f £( 4 )*i + l] + 2^(2)[xj +1 + | ^(3)(xi +1 


+h £< 2 >[vi + l f >(vi + | f > x w)J}) 

+ ±£(3) 4 f x (1 >] +h2 | f J< 3 ) £<«> 4 f x <2) ^ (3) 


f *I< 2 >] + h3 (j 2 £< 2 > £<» f I< 4 > + ^ £« f x< 2 > £<*>] 




} 


X i+1 
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APPENDIX C — Continued 


where 


yi ) 4 

f x (2) = f x (3) A - £j> 

f x (4) A t x (*r +1 ,u 1 * +1 ,t l+1 ) 


(Cl) 


The same equation holds for i = N* - 1, except h is replaced by t£ - jt Q + (N* - l)h]. 
The adjoints for the discrete steepest ascent method are generated by 


AXj = -F?X i+1 


(i = 0, . . N* - 1) 


or 

h = (i + F^x i+ 1 

where 


Ti A i I * * * . . 

F i = lir( x i’ Ui ’ u i + 1 ’ ti ’ ti + 1 


and 


^ r i =x i + l- x i 

as given by the forward Runge-Kutta fourth-order integration. By using the formulas for 
the numerical integration with i = 0, 1, . . N* - 2, 
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APPENDIX C — Continued 

Using this equation yields 



This procedure gives 


\ - {> + h | + 1 f x< 2 > + K <3) + 1 £ (4 > 

+ h 2 | fj(l) fj(2) + | fj(2) £(3) +| fj(3) fj(4) 

+ £< 3 > + T 2 f > f x< 3 > 

+ h4 ^ f xW £< 2 >'x< 3 > £«]}vi 

The same equation holds for i = N* - 1, except h is replaced by t^ - [t 0 + (N* - l)hj. 
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APPENDIX C - Concluded 


Therefore, the only difference between the two methods with respect to the adjoints 
is that is replaced by 



Since x*(t + j - -|j must be approximated from x? and x* + j, the approximations 
might as well be chosen to make f x (2) and ^(3) the same for both methods. 

It has been noted that the replacement for is an approximation for x? +1 

used in the forward Runge-Kutta integration. Likewise, the two approximations (used in 
f x (2) and f x (3)) for x*ft i+ j-^j are used in the forward integration. Thus, one sees 
how the backward integration is being made compatible with the forward integration. 
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APPENDIX D 


COMPATIBLE ADJOINTS FOR PIECEWISE CONSTANT CONTROLS 


In this appendix, the discrete algorithm is simplified for the sample problem by the 
assumption of piecewise constant controls. The form for the difference equations repre- 
senting the state can, in this case, be reduced to 

x i+l ' x i = f i( x i> u i) (i = 0, . . N* - 1) (Dl) 

which is the form discussed at length in reference 3. All changes in the discrete method 
arise from the fact that Hi is now zero. 


Hence, equation (2) becomes 
A(6Xi) ~ F^fiXi + G^fiUi 


(i = 0, . . ., N* - 1) 


where 


9 U 


F ‘ - 9jq( X * ,U i*) 


G i " U i*) 

Equation (3) becomes 

A (x J Oxj) » [x? +1 Fi + Jsxi + A^G^ (1 = 0, 

Equation (4) remains the same; that is, the compatible adjoints of reference 2 are still 
used. Equation (5) becomes 

A(x?’6x i )*xT 1 G i 6u i 

For the sample problem 6x Q is zero; therefore, equation (6) is 
N*-l 

A^j*6 x N * * ^ X S-l G i 6u i 

i=0 

For the sample problem tf is fixed; thus, equations (7) are 




(1 = 0,. . ., N* - 1) 


•> N* - 1) 


= ax N ^N* ,1: N*) 6X N* * ^ > ( X N* ,t N*) " ^ > ( X N* ,t N*) 
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APPENDIX D - Continued 

where di(/ and dfi are omitted since there are no constraints and no stopping condition. 
Equation (8) becomes 

N*-l 

*<*> = 1 x Ii + i G i 6u i <D2) 

i=0 

It is desired to maximize dcp subject to 
N*-l 

(dP)2= ^ SuJ’fiUi (D3) 

i=0 


This procedure is the same as maximizing 


N*-l 


d 4>= I (xT i+1 Gi- M6uT)5u i + /i(dPr 


i=0 

and then solving for the p which will satisfy equation (D3) . Taking the differential yields 

nM 

d(d^- £ xT. +iGi . 2(I6u Tj d(6ui) 
i=0 

This differential will be zero for all d^6ujJ if and only if 


(1 = 0 , 


Substituting the differential into equation (D3) yields 

N*-l 

A 

4 P 

which will be satisfied if 

, 1/2 


<dp)2 =i hi x l 

i*0 ^ 


JL 

2p 


(dP) 2 


<P4> 


which gives for the control corrections 

,lV2 


6u^ = 


(dP) 2 




gFx 


i i+1 


(i = 0, 


N* - 1) 


N* - 1) 
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APPENDIX D — Concluded 


Substituting 6uj into equation (D2) gives 


N*-l 




i=0 


x |,i+l G i 


(dP)2 


)l/2 


G Tx 


<M+1 


“ [W*’ 2 ] 


1/2 


and, as a result, the gradient squared is 
2 



These equations were used in the computer study. 


30 



REFERENCES 


1. Bryson, A. E.; and Denham, W. F.: A Steepest-Ascent Method for Solving Optimum 

Programming Problems. Trans. ASME, Ser. E: J. Appl. Mech., vol. 29, no. 2, 
June 1962, pp. 247-257. 

2. Kelley, Henry J.; and Denham, Walter F.: Modeling and Adjoints for Continuous Sys- 

tems. J. Optimization Theory Appl., vol. 3, no. 3, Mar. 1969, pp. 174-183. 

3. Canon, Michael D. ; Cullum, Clifton D., Jr. ; and Polak, Elijah: Theory of Optimal Con- 

trol and Mathematical Programming. McGraw-Hill Book Co., Inc., c.1970. 

4. Pontryagin, L. S.; Boltyanskii, V. G.; Gamkrelidze, R. V.; and Mishchenko, E. F.: The 

Mathematical Theory of Optimal Processes. Interscience Publ., Inc., c.1962. 

5. Ralston, Anthony: A First Course in Numerical Analysis. McGraw-Hill Book Co., 

Inc., c.1965. 


31 





+2 

10 


>1 


10 


10 ’ 


10 


Id 2 


*©• a- ,«*l 

T3 -o 10 


CM ~ 

I 10 


-o 

o 

W 

o 


-B\ 

10 


10 


10 


10 


10 


10 


B 


-i\ 


O Continuous 1^ (Initial (dP) = 1.0) 
□ Discrete 21 ^ (Initial (dP) 2= 20) 


ga Q 
o 


□ Oo 
□ 


□ 


©□ 


□ 

o 


l°l I I I I I 1 I I I _ L 


J I I 


0 2 4 *6 8 10 12 14 16 18 20 22 24 26 28 30 

Iterations 


Figure 2.- Square of gradient plotted against iteration count (Canon-Cullum-Polak 
algorithm for discrete). Input control, 70°; final time, 10 seconds; computing 
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Figure 3.- Square of gradient plotted against iteration count. Input control, 70°; 
final time, 10 seconds; computing interval, 0.5 second. 
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Figure 4.- Square of gradient plotted against iteration count. Input control, 70°; 
final time, 10 seconds; computing interval, 1 second. 
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Figure 6.- Square of gradient plotted against iteration count. Input control, 80°; 
final time, 10 seconds; computing interval, 0.5 second. 
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Figure 7.- Square of gradient plotted against iteration count. Input control, 60°; 
final time, 10 seconds; computing interval, 0.5 second. 
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Figure 8.- Square of gradient plotted against iteration count. Input control, 70°; 
final time, 20 seconds; computing interval, 5 seconds. 
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